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Abstract 

We  consider  the  use  of  a  multigrid  method  with  central  differencing  to  solve  the 
Navier-Stokes  equations  for  high-speed  flows.  The  time-dependent  form  of  the  equations 
is  integrated  v/ith  a  Runge-Kutta  scheme  accelerated  by  local  time  stepping  and  variable 
coefficient  implicit  residual  smoothing.  Of  particular  importance  are  the  details  of  the 
numerical  dissipation  formulation,  especially  the  switch  between  the  second  and  fourth 
difference  terms.  Solutions  are  given  for  two-dimensional  laminar  flow  over  a  circular 
cylinder  and  a  15  degree  compression  ramp. 
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Introduction 

During  the  1980’s  a  wide  variety  of  numerical  schemes  were  investigated  for  solving 
the  Euler  and  Navier-Stokes  equations.  Multistage  time-stepping  schemes  with  central 
differencing  and  multigrid  acceleration  [1,  2,  3]  were  demonstrated  to  be  quite  effective 
in  computing  subsonic  and  transonic  flows  ■»ver  aerodynamic  components  and  config¬ 
urations.  With  the  recent  resurgence  of  interest  in  high-speed  flight  vehicles,  we  now 
need  to  construct  versatile  algorithms  for  hypersonic  flow.  One  must  keep  in  mind  that 
hypersonic  flows  represent  a  formidable  challenge  for  any  flow  solver.  Ir  particular, 
strong  shock  and  expansion  waves  can  occur  in  the  flow  field,  and  thej  can  interact 
with  each  other  and  with  shear  layers  (i.e.,  boundary  layers,  jets,  wakes).  Such  strong 
nonlinear  behavior  and  interactions  can  easily  cause  divergence  of  any  numerical  inte¬ 
gration  procedure.  This  is  especially  true  during  the  initial  phase  of  a  calculation  with  a 
time-dependent  method.  So  even  the  most  successful  algorithms  of  the  last  decade  may 
require  significant  modifications  to  be  effective  for  hypersonic  flows. 

An  initial  effort  [4]  to  apply  a  central-difference  multigrid  algorithm  to  high-speed 
flows  resulted  in  numerical  difficulties  that  prevented  the  calculation  of  two-dimensional 
flows  (i.e.,  blunt  body  and  wedge  type)  with  a  Mach  number  higher  than  about  7.  In 
order  to  compute  such  flows  a  low  Courant-Friedrichs-Lewy  (CFL)  number  was  required. 
Thus  four  and  five  stage  schemes  were  not  practical,  since  there  is  substantial  deterio¬ 
ration  in  the  high  frequency  damping  of  the  scheme  due  to  the  large  reduction  in  the 
CFL  number.  The  CFL  restriction  reduced  the  potential  of  the  scheme  as  a  viscous  flow 
solver.  More  recently  an  algorithm  utilizing  a  semicoarsening  technique,  a  symmetric 
TVD  formulation,  and  a  three  stage  Runge-Kutta  scheme  [5]  was  proposed  and  used 
to  compute  high  Reynolds  number  (laminar)  Mach  10  flow  over  an  airfoil  at  10  degrees 
angle  of  attack.  A  good  resolution  of  the  bow  shock  wave  and  a  reasonable  convergence 
rate  were  obtained.  The  method  of  semicoarsening  considered  required  a  much  more 
complicated  cycle  strategy  than  that  employed  with  standard  multigrid  methods.  In 
addition,  it  appears  to  be  somewhat  cumbersome  to  implement  in  three  dimensions. 

It  is  our  contention  that  standard  multigrid  techniques  can  be  u.  ^d  in  conjunction 
with  central  differencing  to  compute  hypersonic  flows  effectively.  To  achieve  such  success 
with  these  techniques  one  needs  to  give  appropriate  attention  to  both  the  advection 
and  the  dissipative  processes  of  hyperbolic  multigrid.  The  advection  process  provides 
a  mechanism  by  which  long  wave  disturbances  can  be  rapidly  expelled.  In  a  multigrid 
method  with  a  time-dependent  iterative  procedure,  efficiency  is  in  part  derived  from  the 
larger  time  steps  allowed  on  coarser  meshes.  Hence,  it  is  important  that  the  driving 
scheme  of  the  multigrid  method  use  large  time  steps.  The  dissipative  process  is  essential 
in  smoothing  short  wave  disturbances.  With  this  process  the  multigrid  efficiency  is  based 
on  principles  similar  to  that  for  elliptic  equations. 

For  hypei sonic  flows  one  encounters  an  additional  consideration  regarding  the  damp¬ 
ing  of  the  short  waves.  As  the  Mach  number  increases  the  jumps  across  shocks  become 
larger,  and  it  becomes  more  difficult  to  eliminate  these  high  frequency  oscillations.  Thus 
a  considerable  part  of  the  following  discussion  will  concentrate  on  the  smoothing  al¬ 
gorithm.  The  fundamental  features  of  the  multigrid  process  (i.e..  Full  Approximation 
Storage  scheme,  grid  transfer  operators,  fixed  cycle  strategy)  are  fairly  standard.  Other 
aspects,  such  as  type  of  coarse  grid  correction  scheme  and  procedure  for  smoothing  of 
coarse  grid  corrections,  found  crucial  in  the  present  work  will  be  emphasized.  In  this 
paper  we  consider  a  Runge-Kutta  scheme  [6]  as  the  smoother  for  the  multigrid  method. 
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Central  differences  for  spatial  approximations  are  augmented  by  an  artificial  viscosity 
based  on  TVD  principles  [7].  Several  changes  are  made  to  the  numerical  algorithm  so 
that  a  converged  solution  can  be  obtained  for  high-speed  flows..  We  initially  describe 
the  Runge-Kutta  method  for  the  central-difference  scheme  with  numerical  viscosity.  We 
finally  present  some  examples  to  demonstrate  our  conclusions. 

Basic  Scheme 

The  basic  elements  of  the  scalar  dissipation  model  considered  in  this  paper  were 
first  introduced  by  Jameson,  Schmidt,  and  Turkel  [6]  in  conjunction  with  Runge-Kutta 
explicit  schemes.  The  spatial  discretization  is  based  on  central  differences  with  an  addi¬ 
tional  artificial  viscosity,  This  algorithm  has  been  used  by  many  investigators  to  solve 
the  Euler  equations  ntimerically  for  a  wide  range  of  fluid  dynamic  applications.  The 
same  type  of  spatial  discretization  has  been  applied  to  alternating  direction  implicit 
(ADI)  schemes  [8]  and  LU  factored  implicit  schemes  [9].  In  this  section  the  basic  scheme 
is  briefly  reviewed. 

Consider  the  Euler  equations  in  the  form 


-h /x -b  fl'y  —  0,  (1) 

where  the  four-component  vector  of  conserved  variables 

W  =  [p  pu  pv  pEY  1  (2) 

and  /,p  ar^  the  corresponding  flux  vectors.  The  quantity  p  is  the  density,  u  and  v 
are  the  Cartesian  velocity  components,  and  E  is  the  specific  total  internal  energy.  The 
independent  variables  are  time  t  and  Cartesian  coordinates  (x,j/).  If  (1)  is  transformed 
to  arbitrary  curvilin-^ar  coordinates  ^  =  ^(a:,y)  and  7/  =  rf{x,y),  then  we  obtain 


{J-^W),  +  F^  +  Gr,  =  0,  (3) 

where  J~^  is  the  inverse  transforn  -tion  Jacobian,  and 

E  =  fVr,-  9X„,  G  gx^  -  fy^. 

In  ii  ce" -centered,  finite- volume  method,  (1)  is  integrated  over  an  elemental  volume  in 
the  discretized  computational  domain,  and  J~^  is  identified  as  ti  e  ^-olume  of  the  cell. 
Equation  (3)  can  also  be  written  as 

J-^Wt  +  AW^  +  BWr,  =  0, 

where  A  and  B  are  the  flux  Jacobian  matrices  defined  hy  A  —  dFfdW  and  B  =  dGfdW. 

To  advance  the  scheme  in  time  we  use  a  multistage  scheme.  A  typical  step  of  a 
Runge-Kutta  approximation  to  (3)  is 

W^^')  =  -f  -  AD]  ,  (4) 

where  and  are  spatial  differencing  operators,  and  AD  represents  the  artificial 
dissipation  terms.  The  derivatives  of  the  fluxes  are  approximated  by  central  differences. 
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The  dissipation  terms  are  a  blending  of  second  and  fourth  ditferences.  That  is, 

AD^{DI+DI-D‘-D^)W,  (5) 

where 

=  V,  A(]  W,„  (b) 

DtW  =  V<  AfViA(]  W,j,  (7) 

and  are  the  standard  forward  and  backward  difference  operators,  respectively, 

:\ssoc.iated  with  the  ^  direction.  The  variable  scaling  factor  A  is  chosen  as 


(8) 


where  A^  is  proportional  to  the  spectral  radius  of  the  matrix  A.  The  coefficients  and 
are  ;idapted  to  the  flow  and  arc  defined  as  follows; 


max(i/,-_ij,  Ui+2,j), 

I  Pt+l,j  ~  ^Pt.j  ~b  Pt-i,j  I 


= 


Pi+Ij  +  2pi,j  +  p.-ij 


max 


(9) 

(10) 

(11) 


where  p  is  the  pressure,  and  the  quantities  and  are  constants  to  be  specified. 
The  operators  for  the  t]  direction  are  defined  in  a  similar  manner. 

In  this  paper  we  will  also  consider  a  matrix  form  of  the  dissipation  model  just  de¬ 
scribed.  The  model  of  (5)-(ll)  is  characterized  as  a  scalar  formulation,  since  the  dis¬ 
sipation  for  each  discrete  conservation  equation  is  scaled  by  the  same  eigenvalue.  As 
discussed  in  [7],  a  matrix  form  is  obtained  by  replacing  A  with  a  mauix,  so  that  each 
equation  will  be  scaled  by  its  corresponding  eigenvalue.  That  is,  in  the  ^  direction,  the 
\A\  is  substituted  for  the  eigenvalue  scaling  factor,  A,  in  (6)  and  (7).  For  the  rj  direction, 
(  and  |A|  are  replaced  by  t]  and  |B|,  respectively.  A  convenient  form  for  the  matrix  |/1| 
is  defined  in  the  following  way.  Let 


A  —  Diag  [Ai  A2  A3  A3] 


with 


Ai  =  <7  -f  sJalAal  c,  A2  =  9  -  ^af -f  a|  c,  A3  =  q, 
oi  =  02  =  7  =  OlW  +  0,2V. 

and  c  representing  the  speed  of  sound.  Then, 


i-4|  = 


|A3|/+(M^-|A3| 


7-  1 


E,-\- 


~T~\ — 2 
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H  is  the  total  enthalpy,  and  (j)  =  {v?  +  v^)/2.  Notice  the  special  form  of  (A|,  where  each 
row  of  Ej  is  either  a  scalar  times  the  first  row  or  a  scalar  times  the  second  row  when  the 
first  row  contains  only  zeros.  Due  to  this  special  form  for  any  Ai,  A2,  andAs,  an  arbitrary 
vector  X  can  be  multiplied  by  \A\  very  quickly.  That  is,  one  calculates  (^j+i  *“  “j) 

directly,  rather  than  calculate  and  multiply  a  matrix  times  a  vector.  The  matrix 

|j5|  is  computed  in  the  same  way  as  |/4|  by  simply  replacing  ^  with  t}. 

In  practice  one  cannot  choose  Ai,A2,A3  as  giv  n  above.  Near  stagnation  points  A3 
approaches  zero  while  near  sonic  lines  Ai  or  A2  approach  zero.  A  zero  artificial  viscosity 
would  create  numerical  difficulties.  To  prevent  such  problems,  these  values  are  limited 

as  _ 

iAil  =  max(|Ai|,V;p(A)),  p{A)  =  \q\  +  c^al  +  al,  (13) 

IA2I  =  max(|A2|,K/3(A)),  [Asl  =  maxdAsj,  (14) 

where  the  linear  eigenvalue  A3  can  be  limited  difterently  than  the  nonlinear  eigenvalues. 
The  parameters  and  are  determined  numerically,  and  the  value  used  here  is  0.25. 

The  second-difference  term  in  these  dissipation  models  is  nonlinear.  Its  purpose  is  to 
introduce  an  entropy-like  condition  and  to  suppress  oscillations  in  the  neighborhood  of 
shocks.  This  term  is  small  in  the  smooth  portion  of  the  flow  field.  The  fourth-difference 
dissipation  term  is  basically  linear  and  is  included  to  damp  high-frequency  modes  and 
allow  the  scheme  to  approach  a  steady  state.  Only  this  term  affects  the  linear  stability 
of  the  scheme.  Near  shocks  it  is  reduced  to  zero.  For  high  speed  flows  the  switch  (10)  is 
not  very  good  and  does  not  allow  the  multigrid  to  converge.  Instead  we  consider  a  TVD 
variation  of  the  switch  [7]  given  by 

^ — ,  «(2)  =  1/2.  (15) 

■  -  P.jl  +  \pi,j  -  +  e 

With  this  change  and  the  factor  1/2  in  front  of  the  second-difference  dissipation  term, 
the  scalar  equation  becomes  first-order  upwind  near  shocks.  In  the  case  of  the  original 
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u  we  find  that  v  ~  .05  near  shock  waves  in  transonic  flows.  The  parameter  e  must  be 
chosen  carefully  to  prevent  the  switch  from  being  activated  by  noise.  In  fact  we  found 
it  useful  to  take  an  average  of  the  two  versions  for  v.  Hence,  we  use 


_ |P<+l,i  ^P«ii  ~1~  _ 

(1  -  e)  *  (|p<+i,i  -  Rj\  +  \pi,j  -  Pi~i,j\)  +  e  *  {Pi+i,j  +  2pij  +  Pi-i.j)  ’ 


with  c  =  1/2  a  reasonable  compromise.  We  now  no  longer  have  a  free  parameter  for  the 
second-difference  dissipation. 

Several  other  changes  were  made  to  the  scheme  in  addition  to  the  change  to  a  TVD 
switch.  In  the  original  algorithm  the  artificial  viscosity  for  the  energy  equation  was 
based  on  the  total  enthalpy  rather  than  the  total  internal  energy.  For  high  speed  flows 
we  ba^e  the  artificial  viscosity  on  the  total  internal  energy  so  that  in  each  equation  the 
basic  dependent  variable  is  also  used  in  the  artificial  viscosity.  This  is  more  in  line  with 
upwind  schemes.  This  has  previously  been  used  in  central-difference  schemes  [10].  The 
algorithm  no  longer  preserves  a  constant  total  enthalpy  in  the  steady  state  (as  the  Euler 
equations  do),  but  enthalpy  damping  is  not  useful  for  supersonic  flows.  In  most  cases  the 
difference  between  the  two  approaches  is  small  with  each  approach  having  its  advantages. 
The  original  form  seems  to  give  slightly  sharper  shocks,  while  the  other  one  appears  to 
make  the  scheme  more  robust. 

The  form  of  the  dissipation  model  of  the  basic  or  driving  scheme  is  usually  modified 
for  coarse  grid  problems  in  the  multigrid  process.  A  constant  coefficient  second-difference 
dissipation  is  not  only  less  expensive  computationally  but  also  generally  provides  ade¬ 
quate  smoothing  properties.  For  high  speed  flows  we  find  it  necessary,  as  in  [4],  to  append 
a  nonlinear  dissipation  to  the  usual  one.  Here  this  nonlinear  contribution  depends  on 
the  modified  switching  function  of  (16).  We  also  need  to  increase  the  constant  coefficient 
from  the  standard  value  of  1/16  to  a  value  of  1/4. 

In  order  for  the  scheme  to  be  stable  it  is  necessary  to  restrict  the  time  step.  For 
transonic  flows  it  is  sufficient  to  bcise  this  limitation  on  the  inviscid  terms  except  for 
extremely  fine  meshes.  For  higher  speed  flows  we  found  it  necessary  to  include  a  viscous 
correction  to  the  time  step  restriction  even  for  crude  meshes.  We  shall  thus  develop  a 
sufficient  condition  for  stability  for  the  thin-layer  equations..  Hence,  we  introduce  body 
fitted  coordinates  and  then  ignore  all  second  derivatives  except  for  the  rjr]  derivative. 
One  then  obtains  the  linearized  equation 

Wt  +  Aw^  +  BW^  =  CVF,,,  (17) 

where  again  W  is  the  vector  of  conserved  variables,  and  the  tilde  indicates  that  a  matrix 
is  multiplied  by  the  transformation  Jacobian  J.  The  matrices  of  (17)  are  somewhat 
complicated.  By  transforming  to  nonconservative  variables  the  matrices  are  greatly 
simplified,  and  they  have  the  same  eigenvalues.  Moreover,  W,  A,  B,  and  C  are  replaced 

by 

iy*  =  [p  u  V  pf,  (IS) 

A*  =  MAm-\  B*^.MBM-\  C*  =  MCM-\  (19) 

where  the  matrix  M  is  defined  according  to 
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and  a  is  any  independent  variable.  Abarbanel  and  Gottlieb  [11]  have  shown  that  one 
can  simultaneously  symmetrize  all  these  matrices.  They  symmetrize  the  matrices  with 
the  similarity  transformation  determined  by 

■^00  0 

5  -  ^  °  ®  (21) 

■^p-  0  0  1  0  ’  ''  '' 


■^0  0 


r  0  0  0  1 

\np 

—  0  10  0  ,  . 
-  0  0  10  •  ’ 
-pT-r  0  0  -A-- 

L  n/^1  pc 

Using  this  transformation,  the  triangle  inequality,  and  the  condition  for  symmetric  ma¬ 
trices  that  the  spectral  radius  is  equal  to  the  norm,  one  can  easily  show  that  a  sufficient 
stability  condition  is 


The  quantity  1/A<^  is  bounded  by  the  spectral  radius  of  5“ M*5p  given  by 

and  similarly  the  quantity  1/Atr,  is  bounded  by 

K  =  Hy  -b  u»/*|  -f  c^ril  -b  7/2. 

The  matrix  S~^C*Sp  is  given  by 


■  0  0  0 

A  0  +  VxVy 

3/7  0  TI^T]y  Wy-\rm 

0  0  0 


+  VD 


and  its  spectral  radius  is  given  by 


Kiscous  =  + 

In  general  the  first  term  in  the  maximum  will  be  the  larger.  Hence,  we  can  replace 
^vtscous  in  (23)  by  its  upper  bound.  So  the  actual  time  step  {Atact)  is  determined  as 
follows: 


A<ac<  <  -^4  +  Ap  +  +  vl) 


where  N  is  taken  to  be  the  allowable  CFL  number,  and  the  constant  d  is  4.  In  the 
case  of  steady  flows,  one  can  advance  the  solution  at  each  grid  point  with  the  time  step 
determined  from  this  estimate.  This  type  of  time  stepping  provides  a  preconditioning  of 
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the  matrix  for  the  system  of  diiference  equations.  The  preconditioning  relaxes  stiffness 
due  to  variations  in  local  flow  properties. 

all  flow  calculations  in  this  paper  a  five  stage  Runge-Kutta  scheme  with  a  weighted 
evaluation,  as  detailed  in  [12],  of  the  dissipation  terms  on  the  first,  third,  u'-.J  fifth  stages 
is  used.  As  described  above  the  time  step  is  reduced  in  the  boundary  layer  by  including 
the  viscous  contribution  to  the  time  step.  In  addition,  the  time  step  is  redu':''d  near 
shocks  by  including  a  term  that  depends  on  The  reduction  is  constructed  so  that 
there  is  a  CFL  number  of  1  when  v  =\.  It  serves  to  reduce  the  magnitude  of  the  change 
in  the  solution  near  the  shock  wave,  which  exhibits  strong  nonlinear  behavior. 

Implicit  Residual  Smoothing 

Implicit  residual  smoothing  of  the  residuals  is  used  to  extend  the  stability  range  of 
the  basic  time-stepping  scheme.  For  two-dimensional  problems,  the  residual  smoothing 
can  be  applied  in  the  form 


{I  -  ftV4A()(/  -  =  «!”>,  (26) 

{7rC\ 

where  the  residual  '  is  defined  by 

=  ,  m  =  l,5  (27) 

and  computed  in  the  Runge-Kutta  soage  m,  and  is  the  total  artificial  dissipation 

at  stage  m,  and  is  the  final  residual  at  stage  m  after  the  sequence  of  smoothings 
in  the  ^  and  rj  directions.  The  difference  operators  Cc  and  Co  are  associated  with  the 
convection  and  physical  diffusion  terms.  To  derive  the  maximum  stability  extension  for 
the  hyperbolic  problem,  the  implicit  procedure  is  applied  after  each  stage  of  the  Runge- 
Kutta  scheme.  The  coefficients  and  are  variable  and  functions  of  the  spectral  radii 
and  They  can  be  written  a'  lollows; 


(28) 


where  the  ratio  =  Xr,fX^,  and  the  quantity  N/N*  is  the  ratio  of  the  CFL  number 
of  the  smoothed  scheme  to  that  of  the  basic  explicit  scheme  (usually  having  a  value  of 
2).  In  hypersonic  flow  applications  we  found  it  necessary  for  N*  to  be  3.25,  rather  than 
the  value  of  3.75  used  for  transonic  computations.  From  a  linear  stability  analysis,  the 
scheme  with  these  coefficients  is  stable  for  all  mesh  cell  aspect  ratios  when  the  parameter 
ijj  ^  .125  and  N/N*  is  sufficiently  large.  The  practical  limitation  on  the  Courant  numbci 
is  due  to  the  requirement  for  effective  high  frequency  damping.  For  large  N/N*  the 
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high  frequency  damping  of  the  scheme  vanishes.  The  variable  coefficients  are  functions 
of  the  local  mesh  cell  cispect  ratio,  and  thus  the  smoothing  process  is  not  activated 
in  a  coordinate  direction  where  it  is  not  needed.  This  is  important  for  best  possible 
convergence.  For  further  i.'-cussion  of  implicit  residual  smoothing  see  [13]. 

Multigrid  Method 

As  indicated  earlier  the  sa'ient  features  of  the  multigrid  method  considered  here  are 
fairly  standard.  Moreover,  we  spply  the  Full  Approximation  Storage  (FAS)  scheme  of 
Brandt  [14]  to  define  the  equi-'-lent  fine  grid  problem  on  a  coarse  grid.  Coarser  meshes 
are  obtained  by  eliminating  every  other  mesh  line  in  each  coordinate  direction.  The 
grid  transfer  operators  for  the  solution,  residual,  and  coarse  grid  corrections  are  those 
introduced  by  Jameson  [15].  In  particular,  on  the  auxiliary  meshes,  the  solution  is 
initialized  as 


14?  =  (2S) 

where  the  subscript  refers  to  the  mesh  spacing  value,  the  sum  is  over  the  four  fine  grid 
cells  that  compose  the  2h  grid  cell,  and  is  a  cell  volume.  This  rule  conserves  mass, 
momentum,  and  energy.  On  a  coarse  grid,  a  forcing  function  P  is  added  to  the  governing 
discrete  equations  in  order  to  impose  the  fine  grid  approximation.  After  the  initialization 
of  the  coarse  grid  solution,  this  function  is  computed  as  follows: 


A*  =  r  MWi,)  -  (30) 

where  P/i(W/i)  =  ChWh-  Then,  the  time-stepping  scheme  on  the  (m  -1- 1)^‘  stage  becomes 

=  <>  -  +  4?1-  (31) 

We  can  also  define  a  new  value  W  for  the  residual  as 


R2h  =  R2k{W2H)  +  P2h.  (32) 

collect  this  value,  restrict  the  solution  W2h  to  the  next  coarser  grid,  and  repeat  the 
process.  The  corrections  computed  on  a  coarse  grid  are  transferred  back  to  a  finer  grid 
with  bilinear  interpolation.  In  order  to  execute  the  multigrid  strategy  we  employ  a  fixed 
W-type  cycle.  To  provide  a  well  conditioned  starting  solution  for  the  fine  mesh  a  Full 
Multigrid  (FMG)  method  is  used.  The  FMG  is  analogous  to  grid  sequencing,  except 
multigrid  cycles  are  performed  on  each  coarse  grid. 

Some  of  the  additional  elements  of  the  multigrid  method  are  not  necessarily  standard. 
A  smoothing  of  the  coarse  grid  corrections  being  transferred  to  the  finest  grid  was  found 
to  be  beneficial  in  transonic  computations  [12].  The  smoothing  was  accomplished  with 
the  implicit  residual  .smoothing  mentioned  previously  and  a  constant  coefficient  13  0.1. 

This  smoothing  of  the  residuals  on  the  way  to  finer  meshes  is  crucial  for  the  convergence 
of  the  multigrid  for  hypersonic  flows.  Such  a  proce.ss  acts  to  reduce  high  frequency 
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oscillations  caused  by  the  interpolation.  Hence,  it  becomes  especially  important  near 
strong  shocks,  where  nonphysical  upstream  influence  can  occur.  Another  important 
element  for  high  Mach  number  (M  >  10)  flows  is  the  coarse  grid  correction  scheme. 
That  is,  the  physical  viscous  terms  should  also  be  computed  on  the  coarse  meshes. 

Boundary  Conditions  and  Initialization 

At  a  solid  surface  (wall)  boundary  the  no-slip  condition  is  enforced.  The  wall  pressure 
is  set  to  the  value  at  the  first  interior  solution  point,  and  thus,  a  reduced  normal  mo¬ 
mentum  equation  is  satisfied.  The  wall  temperature  (T^)  is  specified.  In  a  finite-volume 
formulation,  this  amounts  to  treating  the  Cartesian  velocity  components  and  the  tem¬ 
perature  difference  T  —  Tw  scs  antisymmetric  functions  with  respect  to  the  wall.  For  each 
of  the  physical  problems  considered  the  Mach  number  at  the  inflow  boundary  exceeds 
1.0.  Consequently,  the  dependent  variables  are  specified  at  this  boundary  according  to 
the  flow  conditions.  At  any  outflow  boundary,  we  apply  simple  extrapolation  of  the 
components  of  the  solution  vector.  In  general,  for  hypersonic  flows  numerical  difficulties 
are  experienced  at  the  start  of  a  calculation  if  the  discrete  flow  field  is  initialized  with 
free-stream  conditions.  To  avoid  these  difficulties  we  apply  the  following  procedure.  The 
Mach  number  of  the  flow  is  set  to  a  lower  value  (i.e.,  2.0)  than  the  required  one.  In 
addition,  the  wall  temperature  Tyj  is  set  to  the  free-stream  value.  Then  the  Mach  num¬ 
ber  and  Tw  are  gradually  increased  over  a  few  hundred  time  steps  until  the  desired  flow 
conditions  are  obtained.  This  Mach  number  ramping  is  only  done  on  the  coarsest  mesh 
in  the  FMG  sequence. 

Results 

We  consider  two-dimensional  (2-D)  hypersonic  laminar  flow  over  two  different  geome¬ 
tries  in  order  to  evaluate  the  present  multigrid  method.  The  first  geometry  is  a  circular 
cylinder.  For  this  case  the  free-stream  Mach  number  (Mc^)  is  6.5,  and  the  Reynolds 
number  (Reo)  based  on  the  cylinder  diameter  D  (0.82  meters)  is  1.04  x  10®.  The  free- 
stream  temperature  {T<x)  is  202®  Kelvin,  and  the  wall  temperature  is  specified  at  294° 
Kelvin.  This  represents  a  fairly  cold  wall  condition  relative  to  the  temperature  after  the 
normal  portion  of  the  bow  shock.  Computed  surface  pressures  and  heat  transfer  rates 
are  compared  with  the  experimental  data  of  Wieting  [16].  The  second  geometry  is  the 
15  degree  compression  ramp  tested  by  Holden  and  Moselle  [17].  For  this  flow  problem 
the  free-stream  Mach  number  is  14.1,  and  the  Reynolds  number  based  on  a  reference 
length  L  (0.44  meters)  is  1.04  x  10®.  The  length  of  the  flat  plate  preceding  the  ramp  is 
L.  The  Too  is  89°  Kelvin,  and  the  wall  temperature  is  296°  Kelvin.  Surface  distributions 
of  pressure  coefficient  (cp),  skin-friction  coefficent  (c/),  and  heat  transfer  coefficient  (c/,) 
calculated  with  the  present  multigrid  method  are  compared  with  experimental  data  of 
[17].  These  coefficients  have  the  standard  definitions.  In  all  calculations  for  both  cases 
we  assume  that  the  working  fluid  (air)  is  thermally  and  calorically  perfect.  Sutherland’s 
law  is  used  to  determine  the  molecular  viscosity. 

Cylinder  Flow 

The  2-D  cylinder  flow  was  computed  on  a  64  x  64  grid  and  a  128  x  128  grid.  In 
figure  1  the  64  x  64  grid,  which  is  a  proper  subset  of  the  128  x  128  grid,  is  shown. 
For  both  grids  the  circumferential  spacing  is  uniform.-  In  the  normal  direction  on  the 
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centerline,  the  mesh  is  clustered  at  the  surface  with  minimum  spacings  of  approximately 
4  X  10““*  D  and  2  x  10“'*  D  for  the  two  meshes.  At  the  circumferential  angle  (0)  of  -90 
degrees  or  -|-90  degrees,  the  normal  mesh  spacings  are  increased  by  nearly  60  percent 
of  the  centerline  values.  This  was  done  to  accomodate  the  boundary-layer  growth  as 
well  as  the  resolution  of  the  inviscid  flow  region.  As  evident  from  figure  2,  the  normal 
spacing  through  the  shock  region  is  uniform.  Convergence  histories,  which  define  the 
variation  of  the  error  with  multigrid  cycles,  corresponding  to  both  grids  are  displayed 
in  figure  3.  The  error  is  measured  as  an  rms  value  of  the  residual  for  the  continuity 
equation.  In  figure  3  one  observes  three  out  of  the  four  levels  of  refinement  in  the  FMG 
procedure.  The  first  level,  which  requires  only  a  couple  of  CPU  seconds,  is  used  for 
the  initialization  (Mach  number  ramping),  and  thus  it  is  not  shown.  There  are  three 
grids  on  both  the  third  and  the  fourth  levels.  On  the  128  x  128  mesh  the  residual  is 
reduced  nearly  6  orders  of  magnitude  in  300  cycles..  This  requires  about  5  minutes  of 
CPU  time  on  a  Cray  YMP.  It  should  be  emphasized  that  for  engineering  accuracy  (i.e., 
residual  reduced  by  3  orders)  the  finest  mesh  calculation  required  about  2  minutes.  Note 
that  when  engineering  accuracy  is  achieved,  there  is  no  appreciable  improvement  in  the 
viscous  solution  accuracy  by  further  residual  reduction. 

In  figure  4  the  computed  surface  distributions  of  pressure  and  heat  transfer  rate  are 
presented.  There  is  very  good  agreement  between  the  predictions  with  the  128  x  128  grid 
and  the  experimental  data.  For  the  64  x  64  grid  the  scaled  heat  transfer  rate  [QlQrej] 
is  overpredicted  for  0  <  —25  degrees  and  6  >  -f25  degrees.  This  indicates  that  opening 
the  normal  mesh  spacing  adjacent  to  the  surface  produced  a  spacing  too  large  for  the  64 
X  64  grid,  when  only  a  first-order  approximation  is  used  for  the  temperature  derivative. 
The  Mach  number  and  pressure  contours  for  the  two  calculations  are  shown  in  figures  5 
and  6,  respectively.  The  smoothness  of  the  contours  is  evident,  and  the  improved  shock 
resolution  with  mesh  refinement  is  readily  seen.  In  addition,  one  can  notice  that  the 
boundary  layer  is  extremely  thin  for  this  case. 

Compression  Ramp  Flow 

The  2-D  compression  ramp  flow  was  computed  on  grids  consisting  of  56  x  64  (number 
of  streamwise  cells  x  number  of  normal  cells)  and  112  x  128  cells.  Figure  7  depicts  the 
56  X  64  grid.  There  is  stieamwise  clustering  at  the  leading  edge  of  the  flat  plate  and 
at  the  start  of  15  degree  ramp.  The  minimum  spacing  is  approximately  5.8  x  10“^  L. 
Again,  to  resolve  the  boundary  layer,  the  mesh  is  clustered  in  the  normal  direction  near 
the  surface.  The  normal  spacings  for  the  coarse  and  fine  grids  are  about  4.6  x  10““*  L 
and  2.3  x  10““*  L,  respectively.  In  figure  8  the  convergence  histories  for  this  case  with 
the  scalar  and  matrix  forms  of  the  dissipation  model  are  presented..  The  convergence 
rates  with  both  forms  are  quite  good  on  the  56  x  64  grid,  allowing  the  residual  to  be 
decreased  about  6  orders  of  magnitude  in  300  cycles.  The  average  rate  of  reduction  of 
the  residual  with  the  scalar  model  (0.954)  is  slightly  faster.  There  is  some  deterioration 
in  the  rates  with  mesh  refinement.  This  slowdown  with  mesh  refinement  is  also  observed 
for  transonic  computations.  Although  the  average  rates  of  residual  decay  using  the  two 
dissipation  forms  is  essentially  the  same  (0.965),  the  asymptotic  rate  is  faster  using  the 
scalar  model.  The  residual  is  reduced  just  about  5  orders  in  300  multigrid  cycles  with 
the  scalar  formulation.  It  should  be  pointed  out  that  in  the  multigrid  calculation  with 
the  scalar  model  four  grids  were  applied.  Only  three  grids  were  used  in  conjunction 
with  the  matrix  model,  due  to  numerical  difficuHies  caused  by  the  sudden  switch  from 
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free-stream  conditions  to  a  wall  bounded  flow  at  the  inflow  boundary. 

To  provide  a  better  understanding  of  the  computed  results  for  this  compression  ramp 
case,  we  will  first  describe  the  physics  of  the  flow.  Due  to  the  development  of  the 
boundary  layer  on  the  flat  plate,  the  inviscid  flow  is  turned,  and  an  oblique  shock  wave 
is  produced.  Compression  waves  are  formed  by  the  turning  of  the  flow  at  the  start  of  the 
compression  ramp.  These  waves  coalesce  to  form  another  oblique  shock.  The  shock  wave 
emanating  from  the  leading  edge  of  the  plate  intersects  the  compression  ramp  shock.  It 
needs  to  be  emphasized  that  accurate  predictions  of  this  flow  field  depend  strongly  on  a 
good  resolution  of  the  leading  edge  shock.  Also,  resolution  of  the  boundary  layer  on  the 
ramp  is  demanding.  There  is  a  substantial  thinning  of  the  boundary  layer  on  the  ramp 
as  a  consequence  of  the  flow  being  compressed. 

In  figures  9-11  comparisons  are  made  between  the  computed  variations  of  the  pressure, 
skin-friction,  and  heat  transfer  coefficients  and  the  corresponding  experimental  data.  To 
further  aissess  the  shock  capturing  capability  of  the  present  central-difference  scheme, 
results  calculated  with  the  code  developed  by  J.  L.  Thomas,  which  is  based  on  the 
Riemann  solver  of  Roe  and  described  in  [18],  are  also  included  in  these  figures.  The 
computed  distributions  exhibit  excellent  agreement  with  the  data  in  nearly  all  cases. 
With  the  scalar  dissipation  model,  there  are  differences  beween  the  solutions  on  the  56 
x  64  grid  and  the  112  x  128  grid.  The  results  obtained  using  the  matrix  model  for  these 
two  grids  almost  coincide.  Moreover,  the  solution  computed  with  the  matrix  model  on 
a  56  x  64  mesh  is  comparable  to  the  one  calculated  with  the  scalar  formulation  on  a  112 
X  128  mesh.  Figures  12  and  13  show  the  pressure  contours  on  the  ramp  for  each  of  the 
present  computations.  One  can  clearly  see  the  effects  of  dissipation  and  mesh  size  on 
the  leading  edge  shock  and  the  interaction  region  of  the  two  shocks. 

Concluding  Remarks 

A  multigrid  method  with  central  differencing  has  been  successfully  applied  to  the 
solution  of  hypersonic  viscous  flows.  An  explicit  five  stage  Runge-Kutta  scheme  has  been 
used  as  a  smoother  in  solving  the  time-dependent,  thin-layer  Navier-Stokes  equations. 
In  this  paper  considerable  empheisis  has  been  focussed  on  the  dissipative  characteristics 
of  the  driving  scheme  for  the  multigrid  process.  The  presence  of  strong  shocks  has 
required  the  introduction  of  a  switching  function  for  the  numerical  dissipation  based 
on  TVD  principles.  In  addition,  as  a  consequence  of  the  strong  shocks,  a  nonlinear 
coefficient,  which  is  dependent  on  this  switching  function,  has  been  included  in  the 
coarse  grid  dissipation  formulation.  This  nonlinear  coefficient  is  not  needed  for  transonic 
computations.  We  have  also  considered  both  scalar  and  matrix  forms  of  the  dissipation 
model. 

Numerical  solutions  have  been  obtained  for  hypersonic  laminar  flow  over  a  2-D  cylin¬ 
der  and  a  2-D  compression  ramp.  The  agreement  between  predictions  and  experimental 
data  is  quite  good.  Engineering  accuracy  has  been  obtained  rapidly  in  all  computations, 
requiring  about  2  CPU  minutes  on  the  Cray  YMP. 
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Figure  2  Blowup  of  64  x  64  grid  for  2-D  cylinder  flow 


Figure  3  Convergence  histories  for  2-D  cylinder  flow 
computations  (Moo  =  6-5,  ct  =  0.0  deg,  Reo  =  1.04x10®) 
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(b)  Heat  transfer  rate  distribution 


Figure  4  Surface  distributions  of  pressure  and  heat  transfer  rate 
for  2-D  cylinder  flow  (Mo©  =  6.5,  a  =  0.0  deg,  Reo  =  1.04x10® 
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Figure  5  Mach  number  contours,  with  AM  =  0.25,  for  2-D 
cylinder  flow  (Moo  =  6.5,  a  =  0.0  deg,  Rej)  =  1.04  x  10®) 
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Figure  6  Pressure  contours  for  2-D  cylinder  flow 
{Moo  =  6.5,  a  =  0.0  deg,  Rej)  =  1.04x10®) 
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Figure  8  Convergence  histories  for  2-D  compression 
ramp  flow  (Mqo  =  14.1,  Rei  =  1.04x10®,  15  degree  ramp) 
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(a)  Scalar  dissipation 


(b)  Matrix  dissipation 


Figure  9  Pressure  coefficient  distributions  for  2-D  compression 
ramp  flow  (Moo  =  14.1,  Rei  =  1.04  x  10®,  15  degree  ramp) 
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Figure  10  Skin-friction  coefficient  distributions  for  2-D  compression 
ramp  flow  (Moo  =  14.1,  Rei  =  1.04  x  10®,  15  degree  ramp) 
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Figure  1 1  Heat  transfer  coefficient  distributions  for  2-D  compression 
^  ramp  flow  (Moo  =  14.1,  Rei  =  1-04x10^,  15  degree  ramp) 
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(b)  Matrix  dissipation 


Figure  12  Pressure  contours  for  2-D  compression  ramp  flow  on 
56  X  64  grid  (Moo  =  14.1,  Rei  =  1.04x10®,  15  degree  ramp) 
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Figure  13  Pressure  contours  for  2-D  compression  ramp  flow  on 
112  X  128  grid  {Moo  =  14.1,  Rei  =  1.04x10®,  15  degree  ramp) 
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